This PR implements infrastructure for working with matrix sketches, in particular through Gaussian sampling.
The main additions are left_sketch! and right_sketch!, which could be thought of as left_orth with a sketching truncation strategy (although not implemented that way).

Additionally, this brings the CUSOLVER_Randomized algorithm implementation on equal footing by using a Driver to select between the CUSOLVER implementation and the Native one.

Sketching API

• left_sketch[!](A; howmany, kwargs...) -> (Q, B) with Q an m×k isometry and B = Q'·A.
• right_sketch[!](A; howmany, kwargs...) -> (B, Pᴴ) with Pᴴ a k×n co-isometry and B = A·Pᴴ'.
• abstract type SketchingStrategy <: AbstractAlgorithm — supertype for sketching primitives.
• GaussianSketching(howmany; numiter=2, rng=default_rng()) — the one concrete strategy currently shipped.

Sketched SVD

• SketchedAlgorithm{A,S,T,D} — @kwdef wrapper bundling alg (small dense factorization), sketch (SketchingStrategy), trunc (TruncationStrategy), and driver (Driver), to denote decomposing via sketch + decompose + truncate. Mimics TruncatedAlgorithm.
• svd_trunc[!] and svd_trunc_no_error[!] now dispatch on SketchedAlgorithm.
  They allocate (U, S, Vᴴ) sized for the sketch (k = sketch.howmany), check shapes, then call gesvdr!(driver, A, S, U, Vᴴ; sketch, alg, trunc).
• Truncation error: ε = sqrt((‖A‖+‖S‖)(‖A‖−‖S‖)), identical to the old CUSOLVER_Randomized path.
• The CUDA extension now implements gesvdr!(::CUSOLVER, A, S, U, Vᴴ; sketch, trunc, alg=nothing) directly. It currently only accepts sketch::GaussianSketching and trunc::TruncationByOrder (i.e. truncrank) and ignores the inner alg. Requires dedicated overloads for the allocation of the output because of CUSOLVER conventions.
• CUSOLVER_Randomized(; k, p, niters) reimplemented as SketchedAlgorithm(; sketch=GaussianSketching(k+p; numiter=niters+1), truncrank(k); driver=CUSOLVER()). Non-breaking deprecation

Design choices/questions

1. left_sketch! vs left_orth!

I had a hard time deciding between overloading left_orth! directly or introducing a dedicated left_sketch function. For simplicity I kept it separate for now, mostly since this is easier to reason through while fleshing out the design, but there might be a reasonable point to be made to merge the two.

2. Sketching as an algorithm vs Sketching as a truncation strategy

In the svd_trunc implementation, I chose to have both trunc and sketch as a keyword argument, rather than trying to fit this through the same keyword. Again, this is mostly for convenience while I was playing around with this, but we could conceive of a way to merge the two concepts, especially if we decide left_orth!(A; trunc=gaussiansketch(...)) is a reasonable approach.

3. Allocating outputs

One of the things I struggled the most with is the issue of deciding what initialize_output(svd_trunc!, A, ::SketchedAlgorithm) should return.
The main issue is that we already are abusing this a bit for the ::TruncatedAlgorithm codepath, where by convention the USVh we pass in is the one that will be used for the svd_compact! inside, and is not equivalent to the one that is outputted.
For ::SketchedAlgorithm, I therefore chose to pass in the sketched sizes, since that is some workspace that can be reused, and that has similar semantics, but it is definitely true that this is a bit unintuitive.

However, the biggest struggle in changing this is that our automatic differentiation implementations hinge on the fact that everything passes through svd_trunc!(A, out, alg) in the end, so simply not supporting passing in an out argument has quite some implications for the remainder of the code.

This is probably a design choice that we could consider revisiting for a v0.7 version of MatrixAlgebraKit, since it might just be that for these cases where the output size is hard to determine beforehand, it doesn't necessarily make too much sense to allow passing it in as inputs.
The alternative would be to start playing around with "reallocation" strategies, where we shrink/expand provided inputs if they don't have the correct sizes. This does however feel like it might just be a bit too convoluted for what it buys us.

4. Truncation error formula

ε = sqrt((‖A‖ + ‖S‖) * (‖A‖ − ‖S‖))

This is the ‖A‖_F^2 − ‖S‖_F^2 identity, written for numerical stability.
However, I had to significantly lower the tolerances for any such tests, since the numerical accuracy of the sketching and the floating point errors are quite finicky for this.
This is of course the same discussion we had before, where it is not clear to me that svd_trunc really should have the responsibility to provide truncation errors, and svd_trunc_no_error is a way more sensible default mode.
Presumably this is fine as long as it is clearly documented?

5. Gauge fixing depends on the exact truncation path

In the current implementation, the gauge-fixing happens at the moment of decomposing the inner small core tensor, which will then still get unprojected onto a larger space through the sketch.
While it is somewhat straightforward to change this, it might be confusing that there are gaugefix keywords at multiple places, which don't have the same effect.
I'm not sure what the best way forwards is for this, but this might require some more discussion.

TODO

• Documentation page
• Resolve discussions above
• Testing in the context of AD
• Testing in the context of GPU